#Watch out, to reproduct analysis, you have to update the path.
data=read.table("~/Dropbox/Publi_Fusariose/ANALYSIS_REPRO/DATA/bilan_simple_marker.gz" , header=T , sep="," )
#data=data[ order(data$LG , data$Distance) , ]
# random lines of the dataset to develop quicker
#data=data[ sample( seq(1,nrow(data)) , 30000 )  , ]
#data=data[ order(data$LG , data$Distance) , ]

# Genotype matrix
geno = read.table(file="~/Dropbox/Publi_Fusariose/ANALYSIS_REPRO/DATA/fichier_genotypage_QTL.csv", sep = ";" , header = F, na.strings = "-")
geno= as.matrix(geno)
colnames(geno)= geno[1,]
geno= as.data.frame(geno[-1 , ])
# Remove useless individuals
#geno <- geno [ geno[,1] %in% names(data),]

#Watch out, to reproduct analysis, you have to update the path.
pheno=read.table("~/Dropbox/Publi_Fusariose/ANALYSIS_REPRO/DATA/PHENOTYPE/phenotypage_all_fusa_blup.csv" , header=T , sep=";" )
#numeric
pheno[,-1]=apply(pheno[,-1],2,as.numeric)
#check no redondancy in genotype name
table(pheno$geno)[table(pheno$geno)>1]

# put geno name as rowname
rownames(pheno)=pheno$geno
pheno=pheno[,-1]
# delete columns with only NA
which(apply( pheno , 2 , function(x) all(is.na(x)) )==TRUE)

pheno=pheno[ , ! apply( pheno , 2 , function(x) all(is.na(x)) ) ]

library(RColorBrewer)
library(rmdformats)
library(xtable)
library(plotly)
library(FactoMineR)
library(knitr)

my_colors = brewer.pal(8, "Set2") 
my_colors = colorRampPalette(my_colors)(20)

my_carac=levels(data$variable)
nlevels(data$variable)

my_QTL_thres=3.73

# reps
REP1=my_carac[grep("rep1", my_carac )]
REP2=my_carac[grep("rep2", my_carac )]
BLUP=my_carac[grep("blup", my_carac )]

# experiments
CAP13=my_carac[grep("CAP13", my_carac )]
GRI11=my_carac[grep("GRI11", my_carac )]
GRI13=my_carac[grep("GRI13", my_carac )]
GRI15=my_carac[grep("GRI15", my_carac )]
LEC14=my_carac[grep("LEC14", my_carac )]
BEQ11=my_carac[grep("BEQ11", my_carac )]
BAR14=my_carac[grep("BAR14", my_carac )]

# Type de phéno NEPI, NEPIL, PEPIL, NOT, DON
PEPIL=my_carac[grep("PEPIL", my_carac )]
PEPI=my_carac[grep("PEPI", my_carac)]
PEPI=setdiff(PEPI, PEPIL)
NEPIL=my_carac[grep("NEPIL", my_carac )]
NOT=my_carac[grep("NOT", my_carac )]
DON=my_carac[grep("DON", my_carac )]

# AUDPC et derniere notation
AUDPC=my_carac[grep("AUDPC", my_carac )]
LAST=c("CAP13_PEPI550_rep1","CAP13_NEPIL550_rep1","CAP13_PEPIL550_rep1","CAP13_PEPI550_rep2","CAP13_NEPIL550_rep2","CAP13_PEPIL550_rep2","CAP13_PEPI550_blup","CAP13_NEPIL550_blup","CAP13_PEPIL550_blup","GRI11_NOT550_rep1","GRI11_NEPIL550_rep1","GRI11_NOT550_rep2","GRI11_NEPIL550_rep2","GRI13_PEPI550_rep1","GRI13_NEPIL550_rep1","GRI13_PEPIL550_rep1","GRI13_NOT550_rep1","GRI13_PEPI550_rep2","GRI13_NEPIL550_rep2","GRI13_PEPIL550_rep2","GRI13_NOT550_rep2","GRI13_PEPI550_blup","GRI13_NEPIL550_blup","GRI13_PEPIL550_blup","GRI13_NOT550_blup","LEC14_PEPI550_rep1","LEC14_NEPIL550_rep1","LEC14_PEPIL550_rep1","LEC14_PEPI550_rep2","LEC14_NEPIL550_rep2","LEC14_PEPIL550_rep2","LEC14_PEPI550_blup","LEC14_NEPIL550_blup","LEC14_PEPIL550_blup","BEQ11_PEPI450","BEQ11_NOT500","BEQ11_NEPIL500","BAR14_PEPI600_rep1","BAR14_NEPIL600_rep1","BAR14_PEPIL600_rep1","BAR14_PEPI600_rep2","BAR14_NEPIL600_rep2","BAR14_PEPIL600_rep2","BAR14_PEPI600_blup","BAR14_NEPIL600_blup","BAR14_PEPIL600_blup")

# agronomical variable
EPI=my_carac[grep("_EPI" , my_carac )]
FLO=my_carac[grep("FLO" , my_carac )]
HEI=my_carac[grep("HEI" , my_carac )]
VER=my_carac[grep("VER" , my_carac )]

# Resistance set for publication = BLUps when avail. DO we have all the needed traits?
a=c(DON, NOT, PEPI, NEPIL, PEPIL)
b=intersect(a,BLUP)
c=intersect(BEQ11, a)
PUBLI=c(b,c,"LEC14_DON")

# start table
bilan=data.frame( matrix(0,5,7) )
colnames(bilan)=c("CAP13", "GRI11", "GRI13", "GRI15", "LEC14", "BEQ11","BAR14")
rownames(bilan)=c("DON?", "PEPI?", "NOT?", "NEPIL?", "PEPIL?")
all=list(CAP13, GRI11, GRI13, GRI15, LEC14, BEQ11, BAR14)
num1=0
num2=0
for(i in all){
  num1=num1+1
  num2=0
  for(j in list(DON, PEPI, NOT, NEPIL, PEPIL)){
   num2=num2+1
   a=intersect(i, j)
   bilan[ num2 , num1]=length(a)
  }}
kable(xtable(bilan))

#Same function but with interactive graphics!
plot_a_interactiveQTL=function(data, select_chromo, my_variable, LOD_threshold, my_ylim){
  data=data
  select_chromo="4B"
  my_variable=HEI
  
  p=data %>%
    filter(LG==select_chromo & variable %in% my_variable) %>%
    group_by(variable, Distance) %>%
    select(Distance, variable, LOD, R2, group_physique, Posi_physique) %>%
    summarise(
      LOD=max(LOD),
      R2=max(R2, na.rm=T),
      Posi_physique=mean(Posi_physique, na.rm=T),
      group_physique=table(group_physique)[1]
    ) %>%
    mutate(my_text=paste("<br>", "Pos phy: ", group_physique , " | ", Posi_physique , sep="")) %>%
    ggplot(aes(x=Distance, y=LOD, color=variable)) +
      geom_line()
 
# interactive
ggplotly(p)
}

#A function to plot a QTL and visualize the LOD scores along the chromosome. Works for several variable on a specific chromosome.
plot_a_QTL=function(data, select_chromo, variable, LOD_threshold, my_ylim){
  
  # order data
  data=data[ order(data$LG , data$Distance) , ]

  # plot the first variable
  a=data[data$variable==variable[1] & data$LG==select_chromo & !is.na(data$LOD) , ]
  plot(a$LOD ~ a$Distance , type="l" , ylim=c(0,my_ylim) , ylab="LOD score" , xlab="position (cM)" , col=my_colors[1] , lwd=1.6 )
  
  # for every other variable, add a LOD score
  num=1
  for( i in c(2:length(variable))){
    num=num+1
    a=data[data$variable==variable[num] & data$LG==select_chromo & !is.na(data$LOD) , ]
    points(a$LOD ~ a$Distance , type="l" , ylim=c(0,my_ylim) , ylab="LOD score" , xaxt="n" , col=my_colors[num] , lwd=1.6 )
    }
  
  # add lod threshold
  abline(h=my_QTL_thres , col="grey")

  # add legend
  legend("topright" , horiz=F , col=my_colors[c(1:length(variable))] , legend=variable , bty="n" , lty=1 , pt.cex=2, lwd=1.6)
  legend("topleft" , horiz=F , col="white" , legend=select_chromo , bty="n" , lty=1 , pt.cex=0, lwd=1.6)

  }

#A function to summarize all the significant QTL of a list of variables. Gives the max LOD score, the IC, the position etc..
summary_table_for_QTL=function(data, selected_variable, LOD_threshold, size_IC){

  # order data
  data=data[ order(data$LG , data$Distance) , ]
  my_colors=brewer.pal(5,"Paired")
  
  # create empty summary table
  bil=data.frame(matrix(0,1000,10))
  colnames(bil)=c("pop","carac","chromo","LOD_max","position","marker","IC","R2","a","diff")
  
  # loop to study every variable & chromosome
  num_line_bilan=0
    num_chromo=0

    for(chrom in levels(data$LG)){
        num_variable=0
        num_col=0
        num_chromo=num_chromo+1

        for(var in selected_variable){
            num_col=num_col+1
            current_data=data[ data$variable==var & data$LG==chrom  , ]
            current_data=current_data[!is.na(current_data$LOD) , ]
            signif_data=current_data[ current_data$LOD > LOD_threshold  , ]
            
            #Si j'ai des marqueurs significatifs
            if(nrow(signif_data)>0){
                
                #Je récupère les infos de ce QTL
                LOD_max=max(signif_data$LOD)
                mark_max=signif_data$marqueur[signif_data$LOD==LOD_max]
                pos_max=signif_data$Distance[signif_data$LOD==LOD_max]
                r2_max=signif_data$R2[signif_data$LOD==LOD_max]
                a_max=signif_data$a[signif_data$LOD==LOD_max]
          my_diff=ifelse( signif_data$moy.A[signif_data$LOD==LOD_max][1] > signif_data$moy.B[signif_data$LOD==LOD_max][1] , "Dic2>Silur" , "Silur>Dic2" )
            
                #Détermination IC --> je bloque la zone max a 30 cM
                in_IC=current_data[ current_data$LOD > (LOD_max-size_IC) & current_data$Distance > (pos_max-30) & current_data$Distance < (pos_max+30)  , ]
                IC_min=min(round(in_IC$Distance,2))
                IC_max=max(round(in_IC$Distance,2))
                
                #Je remplie le tableau bilan
                num_line_bilan=num_line_bilan+1
                bil[num_line_bilan , 1]="DS"
                bil[num_line_bilan , 2]=var
                bil[num_line_bilan , 3]=chrom
                bil[num_line_bilan , 4]=LOD_max
                bil[num_line_bilan , 5]=pos_max[1]
                bil[num_line_bilan , 6]=as.character(mark_max)[1]
                bil[num_line_bilan , 7]=paste(IC_min, IC_max, sep="-")
                bil[num_line_bilan , 8]=r2_max[1]
                bil[num_line_bilan , 9]=a_max[1]
                bil[num_line_bilan , 10]=my_diff
                
                #close 3 loops
                }}}
    
# Clean and print the Table
bil=bil[bil$pop!=0 , ]
bil=bil[ order(bil$pop , bil$chromo , bil$carac) , ]
return(bil)
#close function
}

#A function to represent confidence interval of QTLs on the genetic map.
show_IC_QTL_on_map=function(data, selected_variable, selected_chromosome, LOD_threshold, size_IC){

  # user can pick every variable
  if(selected_variable[1]=="all"){selected_variable=levels(data$variable)}
  
  # order data
  data=data[ order(data$LG , data$Distance) , ]
  my_colors=brewer.pal(12,"Set3")
  my_colors = colorRampPalette(my_colors)(30)
  var_with_QTL=c()
  
  # loop to study every variable & chromosome
    num_chromo=0

    # prepare background of chart 
    nchro=length(selected_chromosome)
    xchro=c(seq(5,5*nchro,5) )
  par(mar=c(0,5,2,1))
  plot(1,1,col="transparent",bty="n", xaxt="n", yaxt="l" , xlab="", ylab="position (cM)" , xlim=c(3,max(xchro)+11) , ylim=c(320,-5) )
  num=0
  for(chrom in selected_chromosome){
    num=num+1
    A=data[data$variable==selected_variable[1] & data$LG==chrom, ]
    points(rep(xchro[num],nrow(A)) , A$Distance , pch=20 , cex=0.8 , col="grey")
    }

    # IC calculation and plot
  for(chrom in selected_chromosome){

        num_variable=0
        num_chromo=num_chromo+1
        num_col=0

        for(var in selected_variable){
          num_col=num_col+1
            current_data=data[ data$variable==var & data$LG==chrom  , ]
            current_data=current_data[!is.na(current_data$LOD) , ]
            signif_data=current_data[current_data$LOD>LOD_threshold  , ]
            
            #Si j'ai des marqueurs significatifs
            if(nrow(signif_data)>0){
              
              # I add it to the list of variables with QTL
              if(!var%in%var_with_QTL){var_with_QTL=c(var_with_QTL,var)}
                
                # Je récupère les infos de ce QTL
                LOD_max=max(signif_data$LOD)
                mark_max=signif_data$marqueur[signif_data$LOD==LOD_max]
                pos_max=signif_data$Distance[signif_data$LOD==LOD_max]

                # Détermination IC --> je bloque la zone max a 30 cM
                in_IC=current_data[ current_data$LOD > (LOD_max-size_IC) & current_data$Distance > (pos_max-30) & current_data$Distance < (pos_max+30)  , ]
                IC_min=min(in_IC$Distance)
                IC_max=max(in_IC$Distance)
                
                # J'ajoute le trait corresponant a ma variable
                num_variable=num_variable+0.2
                lines(c(xchro[num_chromo]+num_variable,xchro[num_chromo]+num_variable) , c(IC_min,IC_max) , col=my_colors[which(var_with_QTL%in%var)], lwd=6)
                num_variable=abs(num_variable)

        #Close 2 loops and a if
                }}}
    
#Ajout légende et nom de chromosome?
text(xchro , rep(-10,3) , selected_chromosome , col="orange")
#legend("bottomright", 14 , 240 , horiz=F , col=my_colors[1:length(var_with_QTL)] , legend=var_with_QTL , bty="n" , lty=1 , lwd=6 )

#Close function
}

#A function that gives IC features: size, # of markers, # of genes...
find_IC_infos=function(min_IC, max_IC, chrom, to_remove){

  # Get the part of the genetic map concerned by this IC:
  IC_map=data[ which(data$LG==chrom & data$Distance>=min_IC & data$Distance<=max_IC & data$LG==data$group_physique & data$variable==levels(data$variable)[1]), c(1:5)]

  # remove markers if needed.
  if(length(to_remove)>0){
    IC_map=IC_map[-which(IC_map$marqueur%in%to_remove) , ]
  }
  
  # Calculate statistics

  # genet
  deb_IC_gen=min(IC_map$Distance)
  end_IC_gen=max(IC_map$Distance)
  size_IC_gen=end_IC_gen - deb_IC_gen
  size_IC_gen=round(size_IC_gen , 2)
  nb_of_mark=nrow(IC_map)

  # physic
  deb_IC_phy=min(IC_map$Posi_physique)
  end_IC_phy=max(IC_map$Posi_physique)
  size_IC_phy=end_IC_phy - deb_IC_phy
  nb_of_gene=nrow(POS[which(POS$group==chrom & (POS$position*1000000)>=deb_IC_phy & (POS$position*1000000)<=end_IC_phy) , ])

  # return a vec with all info
  vec=c(chrom, min_IC, max_IC, size_IC_gen, nb_of_mark, deb_IC_phy, end_IC_phy, size_IC_phy, nb_of_gene)
  return(vec)
  }

#A function to show all the genes present in a IC. The relatioship between physical and genetic position is also available
show_genes_of_IC=function(min_IC, max_IC, chrom, my_lwd, my_cex, to_remove){
  
  # get data for the IC
  IC_map=data[ which(data$LG==chrom & data$Distance>=min_IC & data$Distance<=max_IC & data$group==data$group_physique & data$variable==levels(data$variable)[1] & data$LG==data$group_physique ) , c(1:5)]

  # remove markers if needed.
  if(length(to_remove)>0){
    IC_map=IC_map[-which(IC_map$marqueur%in%to_remove) , ]
  }

  # physical informations
  #deb_IC_phy=min(IC_map$Posi_physique)
  #end_IC_phy=max(IC_map$Posi_physique)
  my_genes_pos=POS[which(POS$chromo_BW==chrom) , ]

  #plot
  plot(IC_map$Distance ~ as.numeric(IC_map$Posi_physique/1000000),
     ylim=c(min_IC, max_IC),
     xlab="physical pos. of markers (Mb)", ylab="genetic position (cM)" , 
     main=paste("chromome ",chrom ," -- IC : ",min_IC," - ",max_IC," cM",sep="" ),
     col=rgb(0.4,0.3,0.6,0.4) , pch=20 , cex=my_cex
     )
  abline(v=my_genes_pos[,3]/1000000 , col="grey", lwd=my_lwd)

  }

# A function to represent a QTL. It makes a boxplot with 2 categories: alleles A and B
boxplot_a_QTL=function(marker, trait, ...){

  # prepare data: merge expression and alleles:
  nucl= geno[,c(1,which(colnames(geno)==marker)) ]
  my_pheno= data.frame(ind=rownames(pheno), pheno[ , trait])
  AA=merge(my_pheno, nucl, by.x=1, by.y=1, all=T)
  colnames(AA)=c("Genotype", "Phenotype", "Allele")
  AA$Allele[which(AA$Genotype=="dic2")]="A"
  AA$Allele[which(AA$Genotype=="silur")]="B"
  AA=na.omit(AA)
  AA$Genotype
  
  # Highlight Dic2 and Silur
  AA$type=factor(ifelse(AA$Genotype=="dic2" | AA$Genotype=="silur","Highlighted",AA$Allele))
  
  # Make plot
  p=ggplot(AA, aes(x=Allele, y=Phenotype, text=Genotype, fill=type, color=type )) +
    geom_boxplot(alpha=0.5) + 
    geom_jitter() + 
    theme(legend.position="none") +
    ggtitle(trait)

  # interactive
  ggplotly(p)
  }

#A function to calculate interaction between 2 QTLs for a trait
# Function to calculate interaction effect between QTLs:
analyse_inter=function(marker1, marker2, trait){

  # prepare data: merge expression and alleles:
  nucl= geno[ , c(1, which(colnames(geno)%in%c(marker1,marker2)) ) ]
  my_pheno= data.frame(ind=rownames(pheno), pheno[ , trait])
  AA=merge(my_pheno, nucl, by.x=1, by.y=1, all=T)
    AA=na.omit(AA)

    # Prepare
    my_mark=paste(AA[,3],AA[,4],sep="-")
    means <- round(tapply(AA[,2],my_mark,mean,na.rm=T) ,2)
    
    # Calculate interaction
    model=lm(AA[,2] ~ AA[,3] * AA[,4])
    sum=summary(model)
    tot_r2=round(sum$r.squared,3)
    inter_pval=round(sum$coefficients[4,4],10)
    
    # Complete summary files
    res=c(trait, marker1, marker2, means, tot_r2, inter_pval)
    return(res)
}

#A function to boxplot 2 QTLs together
# A function to represent a QTL. It makes a boxplot with 2 categories: alleles A and B
boxplot_two_QTL=function(marker1, marker2, trait, ...){

  # prepare data: merge expression and alleles:
  nucl= geno[ , c(1, which(colnames(geno)%in%c(marker1,marker2)) ) ]
  my_pheno= data.frame(ind=rownames(pheno), pheno[ , trait])
  AA=merge(my_pheno, nucl, by.x=1, by.y=1, all=T)
    AA=na.omit(AA)

    # Prepare
    my_mark=paste(AA[,3],AA[,4],sep="-")
    means <- round(tapply(AA[,2],my_mark,mean,na.rm=T) ,2)

  #plot
    par( mar=c(5,5,2,2))
  boxplot(AA[,2] ~ my_mark , 
      medlwd=0, cex.axis=0.6, cex.col="grey", las=1, main="" ,
      ...,
      ylab=trait,  xlab="Bi-locus genotype" , col=my_colors[c(3,6,8,3)] , xaxt="n", boxwex=0.4) 
  my_labels=c( expression(paste('Dic2'["1B"],"-",'Dic2'["5A"])) , 
               expression(paste('Dic2'["1B"],"-",'Silur'["5A"])) ,  
               expression(paste('Silur'["1B"],"-",'Dic2'["5A"])) , 
               expression(paste('Silur'["1B"],"-",'Silur'["5A"])) )
  axis(labels=my_labels , at=c(1,2,3,4) , side=T)
}

# A function to represent a QTL. It makes a boxplot with 2 categories: alleles A and B
boxplot_two_QTL_interactive=function(marker1, marker2, trait, ...){

  # example to develop if needed
  #marker1="Cluster_16778|Contig1|original@1840"
  #marker2="Cluster_9940|Contig1|complementarySeq@104" 
  #trait="BAR14_PEPIL300_blup"

  # prepare data: merge expression and alleles + add dic2 and Silur
  nucl= geno[ , c(1, which(colnames(geno)%in%c(marker1,marker2)) ) ]
  qq=data.frame(c("dic2","silur"), c("A","B"), c("A","B"))
  colnames(qq)=colnames(nucl)
  nucl=rbind(nucl,qq)
  my_pheno= data.frame(ind=rownames(pheno), pheno=pheno[ , trait])
  AA=merge(my_pheno, nucl, by.x=1, by.y=1, all=T)
    AA=na.omit(AA)

    # Prepare
    AA$my_mark=as.factor(paste(AA[,3],AA[,4],sep="-"))
    means <- round(tapply(AA[,2],AA$my_mark,mean,na.rm=T) ,2)
  AA=AA[order(AA$my_mark),]
  AA$text=paste("ind: ",AA$ind,"<br>","value: ",round(AA$pheno,2),sep="")

  # First we calculate a new x axis with the jiiter added to it. Jitter level must be relative to the number of points per level.
  my_prop=summary(AA$my_mark) / nrow(AA)
  new_x=c()
  for(i in c(1:nlevels(AA$my_mark))){
    myjitter<-jitter( rep(i, summary(AA$my_mark)[i]), amount=my_prop[i]/2)
    new_x=c(new_x, myjitter)
    }

# Plot
plot_ly(AA) %>%
    add_trace( x=as.numeric(AA$my_mark), y = ~pheno, color = ~my_mark, type = "box" ) %>%
    add_trace( x=new_x, y = ~AA$pheno,  
        type = "scatter", 
        mode="markers",  
        marker = list(size = 10, color = ifelse(AA$ind=="dic2","blue",ifelse(AA$ind=="silur","red",rgb(0.1, 0.1, 0.1, 0.3))), line = list(color = rgb(0.1, 0.1, 0.1, 0.8), width = 2)),
        text=~text,
        hoverinfo="text"
    )%>%
    layout(title = "",
           hovermode="closest",
         showlegend=F,
         xaxis = list(title = "" , tickfont = list( color="orange", size=26 ),  tickmode="array", tickvals=c(1,2,3,4) , ticktext=levels(AA$my_mark) , showline=T  ),
         yaxis = list (title = colnames(data)[2] , gridwidth=2, zeroline=F)
         
    )

}